home *** CD-ROM | disk | FTP | other *** search
/ WINMX Assorted Textfiles / Ebooks.tar / Text - Mathematics - Numerical Mathematics and Computing (F).zip / simp.f < prev    next >
Text File  |  2002-06-11  |  3KB  |  109 lines
  1. C
  2. C PAGE 187-190: NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, 1985
  3. C
  4. C FILE: SIMP.FOR
  5. C
  6. C ADAPTIVE SCHEME FOR SIMPSON'S RULE (SIMP,ASMP,PUSH,POP,FCN)
  7. C
  8.       PARAMETER (IST=5,IFS=16)
  9.       REAL STACK(IST,3),FSTACK(IFS,3) 
  10.       EXTERNAL FCN
  11.       DATA EPSI/5.0E-5/, LVMAX/4/     
  12.       A = 0.0 
  13.       B = 8.0*ATAN(1.0)     
  14.       CALL ASMP(FCN,A,B,EPSI,LVMAX,STACK,IST,FSTACK,IFS,SUM,IFLAG)  
  15.       PRINT 3,SUM 
  16.       IF(IFLAG .EQ. 0) THEN 
  17.         PRINT 4   
  18.       ELSE
  19.         PRINT 5   
  20.         DO 2 I=1,IFLAG      
  21.           PRINT 6,FSTACK(I,1),FSTACK(I,2),INT(FSTACK(I,3))
  22.     2   CONTINUE  
  23.       END IF      
  24.    3  FORMAT(//5X,'APPROXIMATE INTEGRAL =',E22.14/)   
  25.    4  FORMAT(10X,'WITH NO BAD SUBINTERVALS')    
  26.    5  FORMAT(10X,'WITH BAD SUBINTERVALS:')      
  27.    6  FORMAT(10X,'[',F10.5,',',F10.5,']',2X,'LEVEL =',I5) 
  28.       STOP
  29.       END 
  30.   
  31.       FUNCTION FCN(X) 
  32.         FCN = COS(2.0*X)/EXP(X) 
  33.         RETURN
  34.       END 
  35.   
  36.       SUBROUTINE PUSH(A,B,LEVEL,STACK,IST,I)    
  37.       REAL STACK(IST,3)     
  38.       IF(I .LT. IST) THEN   
  39.         I = I+1 
  40.         STACK(I,1) = A      
  41.         STACK(I,2) = B
  42.         STACK(I,3) = REAL(LEVEL)      
  43.         RETURN
  44.       ELSE
  45.         STOP 'STACK OVERFLOW IN PUSH' 
  46.       END IF      
  47.       END 
  48.   
  49.       SUBROUTINE POP(A,B,LEVEL,STACK,IST,I)     
  50.       REAL STACK(IST,3)     
  51.       IF(I .GT. 0) THEN 
  52.         A = STACK(I,1)
  53.         B = STACK(I,2)
  54.         LEVEL = INT(STACK(I,3))       
  55.         I = I-1 
  56.         RETURN
  57.       ELSE
  58.         STOP 'STACK UNDERFLOW IN POP' 
  59.       END IF      
  60.       END 
  61.   
  62.       SUBROUTINE ASMP(FCN,A,B,EPSI,LVMAX,STACK,IST,FSTACK,IFS,SUM,IFLAG)      
  63.       REAL STACK(IST,3),FSTACK(IFS,3) 
  64.       COMMON /BLK/ ITOP,LEVEL,LMAX    
  65.       EXTERNAL FCN
  66.       ITOP = 0
  67.       LEVEL = 0   
  68.       IFLAG = 0 
  69.       SUM = 0.0 
  70.       LMAX = 2**LVMAX 
  71.       IF(IST .GE. LVMAX+1 .AND. IFS .GE. LMAX) THEN 
  72.         CALL PUSH(A,B,LEVEL,STACK,IST,ITOP)     
  73.         DO 2 LOOP=1,2*LMAX-1
  74.           IF(ITOP .EQ. 0) RETURN      
  75.           CALL POP(A,B,LEVEL,STACK,IST,ITOP)    
  76.           CALL SIMP(FCN,A,B,EPSI,LVMAX,STACK,IST,FSTACK,IFS,SUM,IFLAG)
  77.     2   CONTINUE  
  78.       ELSE
  79.         PRINT 3   
  80.       END IF      
  81.       RETURN
  82.    3  FORMAT(//5X,'NOT ENOUGH WORKSPACE IN STACK OR FSTACK')
  83.       END 
  84.   
  85.       SUBROUTINE SIMP(FCN,A,B,EPSI,LVMAX,STACK,IST,FSTACK,IFS,SUM,IFLAG)      
  86.       REAL STACK(IST,3),FSTACK(IFS,3),X(0:4),F(0:4) 
  87.       COMMON /BLK/ ITOP,LEVEL,LMAX    
  88.       H = B - A 
  89.       H4 = H/4.0  
  90.       DO 2 I = 0,4
  91.         X(I) = A + REAL(I)*H4 
  92.         F(I) = FCN(X(I))
  93.    2  CONTINUE    
  94.       SUM2 = H*(F(0) + 4.0*F(2) + F(4))/6.0 
  95.       SUM4 = H*(F(0) + 4.0*F(1) + 2.0*F(2) + 4.0*F(3) + F(4))/12.0
  96.       IF(ABS(SUM4 - SUM2) .LE. 15.0*EPSI/REAL(2**LEVEL)) THEN       
  97.         SUM = SUM + SUM4
  98.       ELSE
  99.         IF(LEVEL .LT. LVMAX) THEN     
  100.           LEVEL = LEVEL + 1 
  101.           CALL PUSH(X(2),X(4),LEVEL,STACK,IST,ITOP)       
  102.           CALL PUSH(X(0),X(2),LEVEL,STACK,IST,ITOP)       
  103.         ELSE      
  104.           CALL PUSH(X(0),X(4),LEVEL,FSTACK,IFS,IFLAG)     
  105.         END IF    
  106.       END IF      
  107.       RETURN
  108.       END 
  109.